function [beta2, T, T_J, F, part, OIR] = SURPartiesRobustLagvar_IVGMM(Y, X, Z, mean_a, mean_b, RAND, n, k, l)

% First estimation
P=Z*inv(Z'*Z)*Z';   %Projection matrix
beta=inv(X'*P*X)*X'*P*Y;   % Get parameter estimates
U=Y-X*beta;  % Get residuals

clear P;
% Second Estimation

P2=Z*inv(Z'*(kron(eye(n),ones(2,2)).*(U*U'))*Z)*Z'; %Second step robust projection matrix
beta2=inv(X'*P2*X)*X'*P2*Y; %Second step parameter estimates
u2=Y-X*beta2;  % Second step residuals
U2=(reshape(u2',2,n))'; % Reshape second residuals
alpha=inv([X(:,2) X(:,k+2)]'*[X(:,2) X(:,k+2)])*[X(:,2) X(:,k+2)]'*u2; % Get alpha parameters for error term
ETA=u2-[X(:,2) X(:,k+2)]*alpha; % Get exogenous residuals
clear U P2;

% Compute COVAR matrix
COVAR=inv((X'*Z)*(inv(Z'*(kron(eye(n),ones(2,2)).*(u2*u2'))*Z))*(Z'*X));

% Conduct Wald test on beta parameters
[F, T, T_J] = Wald_GMM(beta2, COVAR, k); 

% Estimate partial effects
bpartial=[beta2(2:k); beta2(k+2:2*k)]; % Create beta vector excluding constant
mean_la=[mean_a;mean_a;mean_a];    %create matrices for partial effects of binary variables
mean_ha=mean_la;
mean_ha(1,4)=1;
mean_ha(2,5)=1;
mean_ha(3,6)=1;
mean_lb=[mean_b;mean_b;mean_b];    
mean_hb=mean_lb;
mean_hb(1,4)=1;
mean_hb(2,5)=1;
mean_hb(3,6)=1;

[part] = partialLagvar_IVGMM(beta2, bpartial, mean_a, mean_b, mean_la, mean_ha, mean_lb, mean_hb, RAND, ETA, alpha, n, k);

% Compute overidentification restrictions test statistic

OIR= (Z'*u2)'*(inv(Z'*(kron(eye(n),ones(2,2)).*(u2*u2'))*Z))*(Z'*u2);
